import numpy as np
import os
import matplotlib.pyplot as plt
import time

def readAccelData(file):
    data = np.loadtxt(file, skiprows=13)
    with open(file, 'r') as csv_file:
        for i in range(10):
            row = csv_file.readline()
        dt = float(row.split()[1])
    return data, dt

file=r'F:\刘中欢\高层所\同门工程\吴梓楠\反应谱计算速度对比\NGA_no_175_H-E12140.txt'
data,dt=readAccelData(file)
# 读取地震波

ndata=len(data) #总点数

T_wave=np.array([dt*i for i in range(ndata)]) #持续时间的各个时间点

a_th = data
T=np.array([0.01*i for i in range(1,1001)])
omega=2*np.pi/T
# funcA可求某频率的单自由度体系在地震作用下的最大加速度,采用newmark-beta法
# omega指单自由度体系的自振频率，para是积分常数，a_th是地震加速度时程，nt时间点的个数（不包括0点）
def funcA(omega, para, a_th, nt):
    # 形成刚度矩阵、质量矩阵和阻尼矩阵，这里是单自由度，假设M=1
    dr = 0.05
    K = omega ** 2
    M = 1
    C = 2 * omega * dr
    P = -a_th
    # 存储各时间点的位移、速度以及加速度
    u = [0]
    v = [0]
    a = [0]
    # 形成等效刚度矩阵K_
    K_ = K + para[0] * M + para[1] * C
    # 递推求解
    for i in range(nt):
        # 计算ti+1时刻的等效荷载
        P_ = P[i] + M * (para[0] * u[-1] + para[2] * v[-1] + para[3] * a[-1]) + C * (
                para[1] * u[-1] + para[4] * v[-1] + para[5] * a[-1])
        # 计算ti+1时刻的位移
        u.append(P_ / K_)
        # 计算ti+1时刻的加速度和速度
        a.append(para[0] * (u[-1] - u[-2]) - para[2] * v[-1] - para[3] * a[-1])
        v.append(v[-1] + para[6] * a[-2] + para[7] * a[-1])
    Aa = []
    for i in range(nt):
        Aa.append(a[i + 1] + a_th[i])
    return np.max(np.abs(u)),np.max(np.abs(Aa))

# 生成某条地震波下的加速度反应谱
# a_th表示经过处理所得地震加速度，omega表示所需的频率范围，funcA可求某频率的单自由度体系在地震作用下的最大加速度
def spec_(a_th, omega, interval, func=funcA):
    # 得到时间点数（不包括0点）
    nt = len(a_th)
    # 选择时间步长t，参数β和γ，并计算积分常数
    t = interval
    gama = 0.5
    beta = 0.25
    para = []
    para.append(1 / beta / t ** 2)
    para.append(gama / beta / t)
    para.append(1 / beta / t)
    para.append(1 / 2 / beta - 1)
    para.append(gama / beta - 1)
    para.append(t / 2 * (gama / beta - 2))
    para.append(t * (1 - gama))
    para.append(gama * t)
    # 计算反应谱
    Sa = []
    for each in omega:
        Sa.append(func(each, para, a_th, nt)[1])
    return (Sa)

start=time.time()
for i in range(1,2):
    Sa=spec_(a_th,omega,dt)
end=time.time()
print(end-start)

plt.plot(T,Sa)
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.title('反应谱') #标题
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
plt.xlabel(r'T/s')
plt.ylabel(r'cm/s^2')
plt.show()